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[i] Phenological dynamics and their related processes strongly constrain land-atmosphere 
interactions, but their relative importance vis-a-vis meteorological forcing within general 
circulation models (GCMs) is still uncertain. Using an off-line land surface model, we 
evaluate leaf area and meteorological controls on gross primary productivity, 
evapotranspiration, transpiration, and runoff at four North American sites, representing 
different vegetation types and background climates. Our results demonstrate that compared 
to meteorological controls, variation in leaf area has a dominant control on gross primary 
productivity, a comparable but smaller influence on transpiration, a weak influence on total 
evapotranspiration, and a negligible impact on runoff. Climate regime and characteristic 
variations in leaf area have important modulating effects on these relative controls, which 
vary depending on the fluxes and timescales of interest. We find that leaf area in energy- 
limited evaporative regimes tends to exhibit greater control on annual gross primary 
productivity than in moisture-limited regimes, except when vegetation exhibits little inter- 
annual variation in leaf area. For transpiration, leaf area control is somewhat less in energy- 
limited regimes and greater in moisture-limited regimes for maximum pentad and annual 
fluxes. These modulating effects of climate and leaf area were less clear for other fluxes and 
at other timescales. Our findings are relevant to land-atmosphere coupling in GCMs, 
especially considering that leaf area variations are a fundamental element of land use and 
land cover change simulations. 
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1. Introduction 

[ 2 ] Numerous climate modeling efforts are focused on 
improving our understanding of the complex interactions 
between the land surface and atmosphere, motivated by 
diverse needs ranging from precipitation generation [Guo 
et al., 2006] to improving our understanding of land use-land 
cover change (LULCC) impacts on biogeochemical cycling 
and climate [e.g., Friedlingstein et al., 2006; Shevliakova 
et al., 2009; Pitman et al., 2011]. To meet these needs, 
dynamic global vegetation models (DGVMs) are increasingly 
being integrated into general circulation models (GCMs) [e.g., 
Foley et al., 1996; Bonan et al., 2003; Zeng et al., 2002; Sitch 
et al., 2003; Krinner et al., 2005; Bonan and Levis, 2006; 
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Shevliakova et al., 2009] with the goal of improving the 
representation of vegetation controls on photosynthesis, 
respiration, transpiration, surface energy fluxes, and carbon 
and nutrient allocation [e.g., Quillet et al., 2010], Yet, a 
recent inter-comparison study on LULCC found significant 
variability in model responses to LULCC, including a lack 
of consistency in how LULCC affects available energy parti- 
tioning between latent and sensible heat [Noblet-Ducoudre 
et al., 2012], These diverse model responses, a partial conse- 
quence of increasing land model complexity, highlight the 
need to better understand the impacts of additional vegetation 
dynamics as they are added to climate models. 

[ 3 ] One of the primary vegetation dynamics categories that 
are poorly understood in a GCM context is vegetation phenology. 
Phenological controls on land-atmosphere interactions in 
GCMs have been studied to clarify the characteristics of 
the connections between phenology and global temperature, 
precipitation, and other climate variables [e.g., Lu and 
Shuttleworth, 2002; Kim and Wang, 2005; Stockli et al., 
2008a], In particular, phenology plays an important role in 
these physical feedback processes through changes in leaf area 
index (LAI) [e.g., Liu et al., 2008; Morisette etal., 2008]. The 
potential influence of phenology on land-atmosphere interac- 
tions is still, however, not well understood, despite the fact that 
the causal link between climate and phenology is conceptually 
straightforward [Alessandri et al., 2007; Morisette et al., 
2008]. Phenological change, for example, can lead to shifts 
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in the surface energy balance and partitioning, which can then 
modify surface temperatures, humidity, and circulation. 

[ 4 ] Presently, GCMs either prescribe phenology (primar- 
ily through specification of leaf area index) [e.g., Friend 
and Kiang, 2005; Lawrence and Chase, 2007] or compute 
phenological variables prognostically [e.g., Sitch et al., 
2003; Krinner et al., 2005; Shevliakova et al., 2009], In 
the former case, phenology is not responsive to the model’s 
surface climatic states (e.g., soil moisture, temperature, and 
solar radiation). In the latter case, phenology does respond 
to climate variations, but the effect of climate biases on 
phenology is then an issue [Kucharik et al., 2006; Morisette 
et al., 2008; Stockli et al., 2008a] making it difficult to pre- 
dict LAI accurately within land models [Liu et al., 2008]. 

[5] Many modeling studies have looked at the influence of 
phenology on land-atmosphere fluxes [e.g., Buermann et al., 
2001; Bonan et al., 2002; Lu and Shuttleworth, 2002; Kang 
et al., 2007; Kim and Wang, 2007], on growing season 
length and carbon fluxes [e.g., White and Nemani, 2003; 
Piao et al., 2007; Richardson et al., 2010; Muraoka et al., 
2010], and on longer-tenn dynamics like tree competition 
[e.g., Kramer et al., 2000]. In fact, many authors have made 
substantial progress in advancing our understanding of the 
relative control of different variables that are important for 
prognostic phenology models. For example, Richardson 
et al. [2007] investigated the relative control of meteorolog- 
ical forcings compared to the biotic response to meteorolog- 
ical variables for a spruce-dominated forest. Alton et al. 
[2007] explored the relative control of biophysical para- 
meters and model processes that govern light propagation 
and canopy photosynthesis for three forest biomes. 

[6] Our interest, however, is on the relative importance 
of phenological and meteorological variables on terrestrial 
carbon and water fluxes in the GCM context. This is impor- 
tant, because many GCMs today still prescribe climatological 
seasonal cycles of vegetation phenology. Given that an 
obvious and imminent development step for these GCMs is 
the inclusion of variable phenology, the impact of such varia- 
tions on the simulation of the surface fluxes that affect simu- 
lated climate must be addressed through careful analysis. 
Our study is a contribution to that analysis. Note that because 
poor simulation of meteorological variables is common in 
GCMs, a determination that meteorological variability is more 
important than phenological variability for fluxes of interest 
may point to meteorological bias correction as a more fruitful 
development path for certain model applications than the 
development of a prognostic phenology routine. 

[7] One key study by Guillevic et al. [2002] made prog- 
ress on this question by exploring the relative controls of 
inter-annual LAI and atmospheric variability on evapotrans- 
piration using the Mosaic land surface model [Koster and 
Suarez, 1992] in both off-line mode and coupled to a 
GCM. The authors found a regional dependence of pheno- 
logical influence on evapotranspiration, where the impacts 
of LAI variability were less both for dense vegetation cover 
(due to a saturation effect) and for regions where vegetation 
experience substantial water stress. They also point out that 
precipitation variability tends to have a stronger impact on 
surface fluxes than LAI and that atmospheric variability is 
the dominate influence in their coupled experiments. 

[8] Here we investigate the relative control of LAI and 
atmospheric variables on terrestrial fluxes of carbon and 


water, building on the work of Guillevic et al. [2002]. Our 
objective is to quantify and compare the strength of LAI 
and meteorological controls on gross primary production, 
evapotranspiration, transpiration, and runoff — all key land 
surface fluxes in GCMs — for a range of climate regimes and 
vegetation. Experiments are designed specifically to assess 
the effects of inter-annual variability in these forcings on 
annual fluxes as well as on maximum values and their timing 
at the sub-annual timescale. We also explore whether the 
contemporaneous nature of LAI and meteorology act together 
to have an impact on carbon and water fluxes. 

[ 9 ] We note that the experiments are conducted with a spe- 
cific land surface model; numerous studies have demonstrated 
substantial variation among the output of land surface models 
due to differences in model parameterizations and forcings 
[e.g., Henderson-Sellers et al., 1995; Dirmeyer et al., 2006; 
Jimenez et al., 2011]. To address this issue, we evaluate model 
performance relative to observations from AmeriFlux sites 
and also highlight important limitations associated with 
this evaluation. 

2. Methods and Model 

[ 10 ] We select four sites in the United States (Figure 1), 
spanning a range of climate regimes and vegetation, to 
explore the relative controls of phenology (through leaf area) 
and atmospheric variables on terrestrial fluxes at spatial 
scales typically associated with global land and climate 
models. The northeast United States (NE) site, located in 
New York state, is covered mainly by deciduous and ever- 
green forest. The southern United States (S) site is mainly 
an agricultural site (primarily including soybean, maize, 
and cotton crops) in Tennessee, while the western United 
States (W) site in western Colorado is a mixture of mainly 
shrubs, woodland, grassland, and evergreen forest. The driest 
location, the far west United States (FW) site in Arizona 
near the southern border with California, is mainly shrubs 
and bare soil. 

[ 11 ] Figure 1 also identifies the locations of three Ameri- 
Flux sites (Harvard Forest, HF; Bondville, B; and Niwot 
Ridge, NR). We select these AmeriFlux sites to validate 
model performance in water and carbon flux simulation, as 
they have observations of latent heat, sensible heat, and 
gross primary productivity fluxes. These AmeriFlux sites 
are used for validation only, because gap-filled meteorolog- 
ical and LAI data are not generally available for as long a 
period (i.e., 17 years) as at our four experimental sites. 

[ 12 ] Figure 2a shows the relationship between annual 
precipitation and temperature, while Figure 2b presents the 
coefficient of variation (CV) of annual precipitation versus 
the standard deviation of annual temperature at our four 
sites. The marker color indicates each site’s evaporative 
fraction (EF), which is the ratio of the site’s annual latent 
heat flux to the sum of its annual latent and sensible heat 
fluxes (computed using the GISS land surface model, as 
described below). EF is a measure of whether evapotranspi- 
ration is soil moisture or energy limited. Low EF values 
(near zero) are indicative of arid and semiarid regions, with 
high-energy inputs and low available soil moisture. High 
EF values (near unity) are typical in humid regions with high 
soil moisture availability, which is generally sufficient to 
satisfy (atmospheric) evaporative demand. 
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Figure 1. Site locations for our experimental analyses plus the three AmeriFlux sites for model valida- 
tion. The experimental sites are 1° by 1° grid cells located in the northeast United States (41.5°N, 
74.5°W), southern United States (35.5°N, 89.5°W), western United States (39.5°N, 106.5°W), and far 
west United States (33.5°N, 114.5°W). The AmeriFlux sites include Harvard Forest (HF; 42.54°N, 
72.17°W), Bondville (B; 40.01°N, 88.29°W), and Niwot Ridge (NR; 40.03°N, 105.55°W). 


[ 13 ] Evapotranspiration is typically energy limited for the 
NE and S sites (EF ss 0.6 for both), meaning that on average 
the soil moisture is high enough to satisfy atmospheric 
demand. Conversely, soil moisture generally controls evapo- 
transpiration at the W site (EF « 0.4), and it is even more 
limiting at the FW site (EF«0.1). With regard to annual 
temperature, the W site is the coolest, while the FW site is 
the wannest. Besides mean amiual climate, we also look at 
the inter-annual variability of temperature and precipitation. 
The standard deviation of annual temperature is similar 
for the four sites («0.6K), but precipitation inter-annual 
variability is somewhat different (Figure 2). The NE and S 
sites have low CV values, while the drier (W and FW) sites 
have larger CV values («0.2 and 0.4, respectively). 

2.1. G1SSLSM 

[ 14 ] All model experiments were conducted with the land 
surface model (LSM) of the NASA Goddard Institute for 
Space Studies (GISS) ModelE, a state-of-the-art atmosphere 
global climate model [ Schmidt et al., 2006; Hansen et al., 
2007], Previous versions of the ModelE have been used to 
investigate soil moisture dynamics [Krakauer et al., 2010], 
the interactions between vegetation and climate [Hasler 
et al., 2009], and the connections between irrigation and 
climate [Puma and Cook, 2010; Cook et al., 2011]. Here 
we model a 1° by 1° grid cell (a resolution that approxi- 
mately corresponds to the new standard resolution of 
ModelE) at each of the four U.S. sites. We use the GISS 
LSM in “off-line mode,” meaning that a data set of meteoro- 
logical forcings (described below) is used to drive the model. 
One advantage to our off-line approach is that we avoid the 
large climate biases typically associated with GCM simula- 
tions; another is that we do not have to deal with uncertainties 
in land-atmosphere feedbacks. 

[ 15 ] The GISS LSM includes both soil and vegetation 
modules that compute water, carbon, and energy dynamics 
of the land surface. The soil module represents the land as 
a 3.5 m soil column, which is discretized into six layers. 


The details of the model’s soil moisture, soil heat, and 
evapotranspiration dynamics are described in Abramopoulos 
et al. [1988], Rosenzweig and Abramopoulos [1997], and 
Aleinov and Schmidt [2006]. The vegetation module of the 
GISS LSM computes photosynthesis using the Farquhar pho- 
tosynthesis model [Farquhar et al., 1980, 1982], Stomatal 
conductance is based on the widely used canopy conductance 
model of Ball et al. [1987]. The approach for representing the 
vertical distribution of canopy radiation is as described in 
Friend and Kiang [2005], which was based on Spitters [1986], 

2.2. AmeriFlux Data 

[ 16 ] The AmeriFlux network, which is part of the global 
FLUXNET network, is a collection of sites with microme- 
teorological towers, which use the eddy covariance tech- 
nique to measure CO 2 , latent heat, and sensible heat fluxes 
[Baldocchi et al., 2001]. FLUXNET sites have been used 
extensively for evaluation, calibration, and parameterization 
of ecosystem and land surface models [e.g., Friend et al., 
2007; Wang et al., 2007; Stockli et al., 2008b; Williams 
et al., 2009]. Here we follow the approach of Stockli et al. 
[2008b] to evaluate the performance of the GISS LSM in 
terms of its predictions of carbon, water, and energy fluxes 
using AmeriFlux observations. 

[ 17 ] We select three AmeriFlux sites (Figure 1) based on 
data availability that approximately cover the range of 
vegetation and climate zone in our experiment (with the 
exception of the arid FW site). The Harvard Forest site is 
characterized as deciduous forest (using the Matthews 
[1984] types) with a sandy loam soil [Urbanski et al., 
2007]. Bondville is cropland, rotated on an amiual basis 
between com and soybean, on silt loam soil [Meyers and 
Hollinger, 2004]. Niwot Ridge, the western AmeriFlux site, 
is an evergreen forest on a clay soil [Monson et al., 2002], 

[is] The atmospheric forcings for these three sites are 
taken from a consistently gap-filled meteorological forcing 
data set developed and generously shared by Stockli et al. 
[2008b], which was based on the original level 2 AmeriFlux 
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Figure 2. Climate characteristics of the four sites based on 
1982-1998 meteorological data from Sheffield et al. [2006]. 
(a) Relationship between annual mean temperature, precipita- 
tion, and evaporative fraction (EF). (b) Relationship between 
coefficient of variation (CV) for annual temperature, CV for 
annual precipitation, and EF. 

data. (For details on the gap-filling methodology, readers are 
referred to Stockli et al. [2008b].) These data are at a temporal 
resolution of 30 min (Bondville, Niwot Ridge) or 60 min 
(Harvard Forest). Beside this meteorological data set, we also 
specify the leaf area index (LAI), because the model’s pheno- 
logical dynamics are prescribed (rather than prognostically 
computed). We use LAI data from the Collection 4 MODIS 
LAI (MOD15A2) product, which is available starting from 
February 2000 [Oak Ridge National Laboratory Distributed 
Active Archive Center, 2011], These 8 day LAI data are 
smoothed using a 5 day moving average and are interpolated 
to match the time step of the atmospheric forcing data. The 
meteorological and LAI data sets allow us to simulate fluxes 
for the years 2001 to 2002 at the Harvard Forest site and the 
years 2002 to 2005 at both Bondville and Niwot Ridge. 

[ 19 ] Our model-simulated fluxes are compared against 
Ameriflux data for latent heat (LE) flux, sensible heat (H) 


flux, and gross primary productivity (GPP) to assess model 
realism in terms of magnitude, timing, and phase of these 
fluxes. Simulated LE and H fluxes are compared against 
level 2 data when the friction velocity (m*) is above 0.2 m, 
which allows us to account for biases in eddy covariance 
measurements during low-turbulence periods [Schmid 
et al., 2003; Stockli et al., 2008b]. The level 2 data set is 
not gap filled, which avoids the imposition of a priori rela- 
tionships on the observations. In the case of GPP, only net 
ecosystem exchange (NEE, which is the difference between 
GPP and ecosystem respiration) data are available in the 
level 2 data set. We therefore use the level 4 data set because 
it includes estimates of GPP. The estimates are obtained 
by applying a “flux-partitioning algorithm” to separate 
NEE into GPP and ecosystem respiration through the use 
of gap-filled NEE (filled by the marginal distribution sampling 
method) and ecosystem respiration (computed using an algo- 
rithm that defines a short-term temperature response of the 
respiration based on NEE) [Reichstein et al., 2005], 

[ 20 ] We account for both random and systematic errors in 
the LE and H observations as done by Stockli et al. [2008b]. 
In particular, we use the empirical work of Richardson et al. 
[2006] to estimate the random uncertainties based on season, 
flux type, and vegetation type. For systematic errors, we 
focus on those associated with the failure in energy balance 
closure. This is done by computing the residual in energy 
balance closure at each site through a regression of observed 
net radiation versus the sum of LE and H fluxes [Wilson 
et al., 2002; Griinwald and Bernhofer, 2007], We then com- 
pute and plot the total (random and systematic) errors for each 
analysis time step (i.e., 10 days) by summing the squares of the 
random and systematic errors and then taking the square root 
of this sum [Stockli et al., 2008b). 

2.3. Experimental Data 

[ 21 ] The atmospheric forcings are from a near-surface 
meteorological data set by Sheffield et al. [2006], which 
were specifically developed for driving land surface models. 
This data set was created by blending reanalysis data with 
observations for the years 1948 to 2008 and is available at 
a 3-hourly, 1° by 1° resolution globally. We select a subset 
of this meteorological data (1982-1998) based on the years 
of LAI data availability. These data are presented in Figure 3, 
where individual years of data are shown by the gray lines 
and the mean monthly values are shown by the blue line. 

[ 22 ] A leaf area index data set is needed for GISS LSM 
simulations, because the model’s phenological dynamics 
are prescribed (rather than prognostically computed). The 
advanced very high resolution radiometers (AVHRR) pro- 
vides the longest continuous global time series of LAI data 
(starting in July 1981) suitable for vegetation studies [ Tucker 
et al., 2005; Ganguly et al., 2008], Although the NASA 
Moderate Resolution Imaging Spectroradiometer (MODIS) 
makes higher quality spectral and angular measurements 
starting from February 2000, it is not available for as many 
years as the AVHRR data. We therefore use the 17 year 
AVHRR total LAI data set (1982 to 1998, excluding 1981 
data) from the International Satellite Land Surface Climatol- 
ogy Project (ISLSCP) Initiative 11 Data Collection [Sietse, 
2010] to take advantage of its longer record. 

[ 23 ] The AVHRR LAI data are biased high when compared 
to MODIS LAI data; consequently, we adjust the AVHRR 
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Figure 3. Monthly (left) temperature and (right) precipitation based on 1982-1998 data for the four sites. 
Mean monthly values are shown by the blue line, while individual year monthly values are shown with the 
gray lines. 
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LAI to remove the high bias, in part because the vegetation 
parameters of the GISS LSM were selected to match observed 
carbon fluxes when using prescribed MODIS LAI data. The 
LAI adjustment is made using a quantile-based mapping tech- 
nique, a method used for many hydrologic and other studies 
[Li et al., 2010], Specifically, the cumulative density function 
(CDF) of the AVHRR data was adjusted to match the CDF 
from the Collection 4 MODIS LAI (MOD15A2) product for 
the years 2001 to 2010. Smooth nonparametric estimates of 
these CDFs were computed using kemal estimation available 
in the MATLAB Statistics Toolbox [Bowman and Azzalini, 
1997]. The nonparametric inverse CDF of the MODIS data 
was then used to generate the adjusted AVHRR LAI data. 
The adjusted AVHRR data are presented in Figure 4 for each 
site; gray lines correspond to individual years of LAI data and 
the blue line corresponds to the mean monthly values. We note 
for the NE site that although LAI during the peak of the grow- 
ing season is not likely to vary dramatically from year to year 
given the characteristics of the forest cover in the region, the 
constant peak LAI is a partial consequence of limitations 
associated with AVHRR measurements. 

[ 24 ] The land cover for each of the four sites is based on the 
1° by 1° land cover map from the ISLSCP II data set [DeFries 
and Hansen, 2010], where land cover is categorized into the 
13 modified International Geosphere Biosphere Program 
classes. This land cover information is then converted to the 


corresponding bare soil and vegetation categories (eight 
vegetation types described by Matthews [1984]) that are used 
by the GISS LSM (Table 1). We note that although the FW 
site is designated to be 19% bare soil, the total bare soil 
fraction is higher, as the shrub category includes interspersed 
bare soil. The shrub and interspersed bare soil are represented 
as a single “tile,” which is the standard approach for the 
GISSLSM and many other global-scale models. 

[ 25 ] As with the vegetation cover and LAI data, soil textures 
(sand, silt, and clay fractions) are also from ISLSCP II 
[Scholes and Brown de Colstoun, 201 1], which contains soil 
parameters on a 1° by 1° grid for two soil depths (0-30 and 


Table 1 . Matthews [1984] Land Cover in Each 1° by 1° Grid Cell 
Based on ISLSCP II Land Cover Data 3 


Land Cover 

NE Site 

S Site 

W Site 

FW Site 

Grassland 

0 

0.02 

0.18 

0.01 

Shrub 

0.04 

0.18 

0.18 

0.74 

Woodland 

0.02 

0.10 

0.33 

0 

Deciduous forest 

0.57 

0.06 

0.01 

0 

Evergreen forest 

0.34 

0.07 

0.23 

0 

Crops 

0 

0.56 

0.08 

0.06 

Bare Soil 

0.03 

0.02 

0 

0.19 


a See DeFries and Hansen [2010], 
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Figure 4. Monthly LAI from the AVHRR global data set (adjusted using MODIS data) for 1982-1998 
data at the four sites: (a) northeast U.S. site, (b) southern U.S. site, (c) western U.S. site, and (d) far west U. 
S. site. Mean monthly values are shown by the blue line, while individual year monthly values are shown 
with the gray lines. 


19 


PUMA ET AL.: CONTROLS ON LAND- ATMOSPHERE FLUXES 


0-150 cm). The soil textures for the four sites are shown 
in Table 2. Other land-related parameters are as described 
in Abramopoulos et al. [1988] and Rosenzweig and 
Abramopoulos [1997]. 

3. Experimental Design 

[ 26 ] The experimental setup, presented schematically in 
Figure 5, is designed to assess the relative control of both 
leaf area variations and meteorological variations on varia- 
tions in terrestrial carbon and water fluxes. For each site, 
we run individual yearlong simulations for every combina- 
tion of the yearly LAI and atmospheric data, which results 
in 289 realizations (i.e., 17 LAI years x 17 meteorological 
years). All experiments are initialized using the climatological 


Table 2. Soil Textures Fractions for Two Soil Depths (0-30 and 
0-150 cm) at the Four Sites Based on ISLSCP II Soil Data 3 


Soil Texture 

NE Site 

S Site 

W Site 

FW Site 

Sand fraction (0-30 cm) 

0.36 

0.40 

0.42 

0.57 

Silt fraction (0-30 cm) 

0.33 

0.36 

0.36 

0.26 

Clay fraction (0-30 cm) 

0.31 

0.24 

0.22 

0.17 

Sand fraction (0-150 cm) 

0.42 

0.40 

0.40 

0.51 

Silt fraction (0-150 cm) 

0.30 

0.33 

0.31 

0.28 

Clay fraction (0-150 cm) 

0.28 

0.27 

0.29 

0.21 


a See Scholes and Brown de Colstoun [2011], The 0-30 cm soil textures 
are assigned to the top two GISS LSM soil layers (thicknesses of 0.1 and 
0. 1 7 m, respectively), and the 0-150 cm data are assigned to the bottom four 
layers (thicknesses of 0.30, 0.51, 0.89, and 1.53, respectively.) 


soil moisture and soil heat values for January based on a 
1 7 year model run with each site’s actual LAI and meteorolog- 
ical data from 1982 to 1998. From these simulations, we com- 
pute two sets of average values to explore the control of LAI 
variability (LAI_v) and meteorological variability (MET_v) 
on carbon and water fluxes. Each member of the first set 
of average values is the average across simulations run with 
a single LAI year but different meteorological years, so that 
the effect of the meteorological variability is averaged out. 
The distribution of these average values ( A| a j| . . . A'lai 1 7 
shown in the blue boxes of Figure 5) reflects the impact 
of varying LAI on carbon and water fluxes. For a given 
member of the second set, the average flux values are 
computed across simulations that use the same meteorolog- 
ical year but different LAI years, so that the effect of LAI 
variability is averaged out. The distribution of this second 
set of average values (X met i . . .W met i 7 presented in the red 
boxes of Figure 5) reflects the impact of varying meteoro- 
logical forcings on terrestrial fluxes. 

[ 27 ] In addition to these two sets of average values, the 
shaded boxes along the diagonal have a special significance, 
because the LAI and meteorological data set are from the 
same year. Therefore, the terrestrial fluxes corresponding 
to these shaded boxes represent the fluxes due to the actual 
historical LAI and meteorological forcings at each site. We 
analyze the distribution of these “actual” fluxes (relative to 
the fluxes from all of the realizations) to investigate whether 
the contemporaneous nature of extreme LAI values and 
meteorological conditions can exacerbate the extremes of 
carbon and water fluxes. 
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Figure 5. Schematic diagram of the experimental setup to assess the relative control of leaf area index 
and meteorological variables on terrestrial carbon and water fluxes for the four experimental sites. 
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4. Results 

4.1. Model Validation 

[ 28 ] The model simulates GPP fluxes well in terms of 
timing; this is especially true at Harvard Forest (Figure 6). 
The additional model skill in terms of GPP fluxes at Harvard 
Forest is partially related to the original vegetation model 
development, which involved selection of biophysical para- 
meters that are both consistent with the site’s vegetation and 
allow accurate simulation of GPP fluxes. At Bondville, the 
small springtime increases in GPP missed by the model are 
not reflected by corresponding increases in MODIS LAI 
data (not shown). Similarly, the delayed springtime green- 
up and early senescence at Niwot Ridge is consistent with 
the MODIS LAI timing. In terms of magnitude, the GISS 
LSM matches the Harvard Forest fluxes very well but under- 
estimates the peak growing season LAI at the other two sites. 
However, full evaluation and attribution of this underestima- 
tion is outside the scope of this paper, as it requires detailed 


knowledge of the uncertainties associated with GPP estimates. 
These uncertainties can be substantial due to both gap-filling 
issues and the assumptions associated with ecosystem respira- 
tion modeling [Falge et al., 2002; Hagen et al., 2006]. 

[ 29 ] Overall, the LE and H fluxes simulated by the GISS 
LSM are in close agreement with observations, as shown 
in Figures 7 and 8. At Harvard Forest, simulated LE fluxes 
are within the range indicated by the uncertainty bars. The 
model is not as successful with the observed seasonality of 
H fluxes (missing the springtime peak), which is a deficiency 
also found in deciduous forest simulations (Morgan Monroe 
State Forest) with the Community Land Model version 3.5 
(CLM3.5) [Stockli et al, 2008b], 

[ 30 ] The GISS LSM is also quite skillful at simulating LE 
flues at Bondville, where the peak and timing of the simulated 
fluxes are in remarkable agreement with observations. The 
model performs well in terms of simulating H seasonality, 
although the magnitude of these fluxes tend to be lower than 
observations. At Niwot Ridge, the seasonality of modeled 


Harvard Forest 



Year 


Bondville 



Figure 6. Comparison of model-predicted and level 4 observed GPP fluxes (10 day averages) at the three 
AmeriFlux sites. Model data are represented by a line with squares; observations are represented by stars. 
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Figure 7. Comparison of model-predicted and level 2 observed latent heat fluxes (10 day averages) at the 
three AmeriFlux sites. Error bars represent the estimated (random and systematic) uncertainties associated 
with the observations. Model data are represented by a line with squares; observations are represented by 
stars. Note that for Harvard Forest, net radiation observations were not available for 2002, so these uncer- 
tainties were not estimated. 


LE fluxes is in agreement with observations, but simulated LE 
fluxes are lower than observations. H fluxes are simulated 
quite well in tenns of timing and magnitude. 

[ 31 ] Part of the discrepancies in flux magnitudes can be 
attributed to the failure in energy balance closure, which 
we estimate using linear regression coefficients [e.g., Wilson 
et al., 2002]. The energy balance closure problem is most 
significant at both Bondville (slope: 0.64; intercept: 19 W/m 2 ) 
and Niwot Ridge (slope: 0.73; intercept: 13 W/m 2 ) and less 
problematic at Harvard Forest (slope: 0.83; intercept: 
—3.9 W/m 2 ); ideal energy balance closure would have a slope 
of 1 and intercept of 0 [Wilson et al., 2002]. Deficiencies in 
model processes are also possible. For example, the wintertime 
LE underestimate suggests possible deficiencies in the snow 
evaporation/sublimation algorithm [Lynch-Stieglitz, 1994], 
given the persistent snow cover at this site during winter 
[Monson et al., 2005]. Growing season discrepancies may 


vary due to vegetation parameters, such as the water stress 
factor, that are not optimized for the vegetation at a particu- 
lar AmeriFlux site. 

[ 32 ] The main goal associated with this AmeriFlux evalua- 
tion is to determine whether the GISS LSM is capable of real- 
istically simulating LE, H, and GPP fluxes. While there is 
room for improvement (as with all state-of-the-art models), 
the results here (Figures 6-8) demonstrate that in general, 
the GISS LSM skillfully simulates these dynamics in terms 
of timing and magnitude. We therefore have confidence in 
applying the model for the experiments described by Figure 5. 

4.2. Experiments 

[ 33 ] The simulation results are presented for GPP, evapo- 
transpiration (ET), and transpiration in Figures 9—11, respec- 
tively. In each, we present a four-paneled figure with box 
plots of the annual fluxes, the timing of the 50% cumulative 
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Figure 8. Comparison of model-predicted and level 2 observed sensible heat fluxes (10 day averages) at 
the three AmeriFlux sites. Error bars represent the estimated (random and systematic) uncertainties asso- 
ciated with the observations. Model data are represented by a line with squares; observations are repre- 
sented by stars. Note that for Harvard Forest, net radiation observations were not available for 2002, so 
these uncertainties were not estimated. 


annual flux, the value of the maximum pentad (5 day) flux, 
and the timing of this maximum pentad flux. As with stan- 
dard box plots, the median is the central line and the 25th 
and 75th percentiles are the box edges. The whiskers extend 
to the lowest datum within 1.5 times the inter-quartile range 
(IQR) of the lower quartile and to the highest datum within 
1.5 times the IQR of the upper quartile. The remaining data 
outside this range are plotted separately with markers. 
Model simulations with historical LAI and meteorological 
forcings are presented at the end of this section. Note that 
for all sites, time is measured from 1 January. Alternatively, 
we could present timing results based on the start of the 
water year (i.e., 1 October) at the western (W and FW) sites, 
but because the results turn out to be the same, we do not 
distinguish starting dates among the sites for clarity. 

[ 34 ] Figure 9 presents the box plots for GPP at the four 
sites. We first note that vegetation productivity, as measured 


by total annual GPP, decreases moving from east to west 
(Figure 9a). At the sites where energy limits evapotranspira- 
tion (NE and S sites), differences in annual GPP are due 
primarily to the vegetation cover (i.e., forest versus crops). 
Reduced productivity by vegetation at the W and FW sites 
is due to the joint impacts of reduced vegetation cover and 
water availability. At all sites, LAI control (LAl_v) on the 
distribution of annual GPP is clearly larger than meteorolog- 
ical control (MET_v), as indicated by the broader range of 
values seen for the LAl_v averages. Strong LAI control 
at the S site stands out, a consequence of the site’s rela- 
tively large inter-annual variability in maximum growing 
season LAI (Figure 4b), which is related to year-to-year 
differences in crop distribution and production of soybean, 
maize, and cotton. 

[ 35 ] We next consider the relative controls of LAI and 
meteorology on the temporal distribution of fluxes within a 
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Figure 9. Box plots of gross primary productivity for the LAI and meteorological experiments: (a) 
annual values, (b) timing of the cumulative 50% of the annual flux, (c) maximum pentad flux, and (d) timing 
of the maximum (pentad) flux. 


year using the timing of the cumulative 50% of the annual 
flux. The variability in this timing is relatively modest (i.e., 
1QR S3 1 pentad) except at the F W site (IQR ~ 7 pentads 
for the LAI_v case). Why is the variability in the timing so 
large at the FW site (Figure 4)? First, it is clear from Figure 9 
that the magnitude of GPP at the FW site is relatively 
small. Also, vegetation productivity at the site is intermittent, 
typically linked to periods of sufficient water availability. 
These intermittent periods of productivity are large enough 
to account for a substantial fraction of the annual total GPP. 
The LAI of the NE site shows variation in timing of spring 
green-up and autumn senescence (Figure 4a); this variation 
results in fairly strong control of LAI on GPP at the site 
(Figure 9a). At the S and W sites, substantial LAI control 
on GPP occurs due to variation not only in green-up and 
senescence timing but also in maximum LAI (Figures 4b 
and 4c, respectively). 

[36] Moving beyond the annual timescale, we next look at 
the maximum pentad (5 day) GPP fluxes (Figure 9c) and the 
timing of these pentad fluxes (Figure 9d). The magnitude of 
these maximum pentad fluxes for the LAI_v decreases 
across sites moving from east to west, in analogy to the 
results for annual GPP. At the NE site, the maximum pentad 
flux is essentially a fixed value, a direct consequence of the 
fixed peak growing season LAI (Figure 4a). On the other 
hand, the maximum LAI at the S site does display substan- 
tial inter-annual variability (Figure 4b), which corresponds 
to significant LAI control on the maximum pentad GPP. 


[ 37 ] The timing of the pentad GPP fluxes is presented in 
Figure 9d to show LAI and meteorological controls on peak 
growing season productivity. Even though the NE and S sites 
are in energy-controlled evaporative regimes, the relatively 
short timescale associated with the maximum pentad flux 
(compared to the 50% cumulative flux) allows meteorological 
forcings to have a greater influence. As expected, LAI at the 
NE site has minimal control on maximum pentad flux timing, 
while LAI control at the S site is only slightly larger than 
the meteorological control. Interestingly, the meteorological 
control at the W site is dominant, while the influences of 
LAI and meteorological forcings on timing of the pentad 
GPP fluxes are comparable at the FW site. 

[38] The most notable difference between GPP and ET is 
that meteorological variability exerts a larger control than 
LAI variability for all quantities in Figure 10. At the NE site, 
the lack of LAI control is a partial manifestation of LAI 
seasonality, where the maximum growing season LAI is 
essentially constant. At the S and W sites, LAI variability 
is not inconsequential, although meteorological effects are 
still larger for each quantity in Figure 10. The strongest 
relative ET control of meteorology occurs at the FW site, 
where bare soil coverage and the associated soil evaporation 
are substantial (19% bare soil). Also, the site’s small LAI 
values diminish any phenological influences on ET. 

[ 39 ] LAI has the greatest relative control on annual ET at 
the S site (Figure 10a), which is linked to the relatively large 
year-to-year differences in LAI. LAI also influences ET 
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Figure 10. Box plots of evapotranspiration for the LAI and meteorological experiments: (a) annual 
values, (b) timing of the cumulative 50% of the annual flux, (c) maximum pentad flux, and (d) timing 
of the maximum (pentad) flux. 


timing, where its strongest relative control on the timing of 
its 50% cumulative flux (Figure 10b) is at the S site and 
slightly less at the NE site. LAI exerts an important control 
on the cumulative flux timing at the W and FW sites 
(although meteorological forcings are still more influential). 
LAI control on the timing of its 50% cumulative flux 
(Figure 10b) is again strongest at the S site and slightly less 
at the NE site. 

[ 40 ] The LAI influence on the maximum pentad flux and its 
timing (Figures 10c and lOd) does play a role at the S site, 
which has an energy-limited evaporative regime and signifi- 
cant inter-annual LAI variability. The NE site has no inter- 
annual variability in peak growing season LAI, while soil 
moisture controls at the FW site are extremely strong. There- 
fore, meteorological controls on the maximum pentad fluxes 
and their timing are dominant with LAI having a very small 
influence. The influence of LAI at the semiarid W site is 
slightly greater, because it is slightly wetter and has greater 
LAI inter-annual variability. 

[ 41 ] Figure 1 1 separates out the impacts on transpiration, 
for which we find that neither LAI nor meteorological 
forcings have a consistently larger influence on all the plot- 
ted quantities. The relative control of LAI on annual transpi- 
ration becomes larger moving from the NE site to the FW 
site (Figure 11a) as does the maximum pentad transpiration 
(Figure 11c). The temporal distribution of transpiration 
within a year (Figure 1 lb) also shows comparable influences 


of LAI and atmospheric forcings, where meteorological 
controls are somewhat stronger, except at the FW site. Here 
the LAI_v experiment has a large IQR for LAI v (6.2 pentads), 
which is interesting in the context of GPP (Figure 9b) and ET 
timing (Figure 10b). That is, while the importance of LAI is 
similar for transpiration and GPP timing, it is unimportant 
for ET timing. This result is explained by the fact that the 
average contribution of transpiration to total ET at this site is 
only 14%. As for timing of the maximum transpiration flux, 
meteorological controls are most important at the energy- 
limited sites, whereas LAI controls are more substantial at 
the FW site (Figure lid). 

[ 42 ] Surface and subsurface runoff are also relevant sur- 
face water fluxes. LAI variability potentially impacts surface 
runoff by modifying transpiration and canopy interception 
of precipitation. However, a box plot (not shown) of LAI 
and meteorological controls on the surface runoff fluxes 
and their timing reveals that LAI control is negligible 
relative to meteorological controls at all four sites. The same 
is true for subsurface runoff. 

[ 43 ] Guillevic et al. [2002] point out that in real-world, 
moisture-limited regimes, periods of lower precipitation 
are often associated with periods of lower vegetation 
density and vice versa. To assess whether our model repro- 
duces this expected behavior, we compare in Figure 12 the 
standard deviation of GPP, ET, and transpiration for the 
17 realizations with the actual historical combinations of 
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Figure 11. Box plots of transpiration for the LAI and meteorological experiments: (a) annual values, 
(b) timing of the cumulative 50% of the annual flux, (c) maximum pentad flux, and (d) timing of the 
maximum (pentad) flux. 


LAI and meteorological forcings ((factual) t0 that obtained 
for the entire set of 289 realizations (u^aii)- If °A,actuai 
exceeds g x ,m, then it provides evidence that correlations 
in nature between LAI and atmospheric variables have 
a compounding effect on the extremes of carbon and 
water fluxes. 

[ 44 ] Figure 12 (top) presents the standard deviation com- 
parison for annual GPP, which shows that OGPP.actuai 
exceeds tfGPP.aii at all sites. In middle panel, we find that 
(7 ET.actuai exceeds Get, ail at three sites, with the W site being 
the exception. The bottom panel shows transpiration 
results. For reference, the percentages of transpiration in 
annual ET at the NE, S, W, and FW sites are approximately 
48%, 42%, 42%, and 25%, respectively. We find that the 
standard deviation for the actual forcings is larger than 
the whole-realization standard deviation except for ET at 
the W site and transpiration at the S site, although the 
differences are not statistically significant based on a two- 
sided Student’s t test. (This is not surprising given the small 
number of actual forcing years.) Then, for the most part, 
our results support the idea that correlations in LAI and 
meteorological anomalies help exacerbate extremes in 
surface fluxes. This is particularly indicated for the NE 
and FW sites, for which an analysis of the histograms of 
the fluxes (not shown) reveal a shift in the distribution of 
the fluxes toward the extremes. 


5. Discussion and Conclusions 

[45] The simulation of vegetation in climate models has been 
the focus of considerable model development and research 
efforts. What remains unclear is the importance of LAI 
variability relative to that of atmospheric variables for 
different climate regimes and vegetation. As modelers make 
decisions with regard to model complexity for representing 
land and vegetation dynamics, a clearer understanding of 
how LAI impacts land-atmosphere interactions is essential. 
We find that the relative importance of LAI phenology 
varies substantially depending on the flux of interest, the 
characteristics of LAI variability, and the climate regime. 
Our analyses indicate that relative to meteorological variables, 
LAI has a dominant control on gross primary productivity, a 
comparable (but smaller) influence on transpiration, a 
weaker influence on evapotranspiration, and a negligible 
impact on runoff. 

[46] GCMs use an energy balance approach to compute 
evapotranspiration and sensible heat fluxes; in this approach, 
surface-to-air temperature and humidity differences are key 
driving forces controlling these fluxes [e.g., Cornwell and 
Harvey, 2007]. In comparison, GPP fluxes are modeled 
through a dependence on relative humidity and leaf surface 
CO 2 concentrations as well as through environmental stress 
factors that affect the vegetation’s photosynthetic capacity. 
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Figure 12. Standard deviation of actual LAI and meteoro- 
logical forcings (o^actuai) compared with the standard devia- 
tion of all simulations (Vr Y , aM ) for GPP, ET, and transpiration. 

Although GPP and the transpiration component of ET are 
coupled through a common canopy conductance term, meteo- 
rological controls on these fluxes are not equivalent. We see, 
for example, that at the NE site, LAI has a dominant control 
on annual GPP, while meteorology has the dominant control 
on annual transpiration. This difference is a consequence of 
the complex interplay of meteorology, LAI, and vegetation 
biophysical properties at the site. 

[ 47 ] The modulating effects of climate regime and LAI 
variations on the relative control of LAI and meteorology 
depend on the flux and timescale of interest. In terms of 
climate effects, evaporative regime (as estimated by evapo- 
rative fraction) is useful as a basic indicator of LAI control 
on annual and maximum pentad GPP; LAI at energy-limited 
sites tends to exhibit greater control than at the moisture- 
limited sites. Other characteristics of LAI (e.g., seasonality 
and inter-annual variability), however, also appear to be 
important modulators of LAI control on surface fluxes. This 


is clearly demonstrated at the NE site, where the maximum 
growing season LAI is constant from year to year, thereby 
reducing LAI control on GPP (Figures 9a and 9c). Evapora- 
tive regime is also valuable for assessing the relative LAI 
control on annual and maximum pentad transpiration, where 
relative LAI control is somewhat less at energy-limited sites 
and greater at moisture-limited sites (Figures 11a and 11c). 
For total ET, LAI has much less impact on variation in ET 
than does meteorology. At shorter timescales that are impor- 
tant for the timing of fluxes, the modulating effects of climate 
regime and LAI are unclear. One exception is at the FW site, 
where precipitation variability is very large and, consequently, 
meteorological controls dominate ET timing. 

[ 48 ] The strong control of LAI on gross primary productiv- 
ity points to the need to have accurate representation of LAI 
variability within climate models. Climate modelers focused 
on the carbon cycle therefore need to select approaches for 
representing LAI phenology, whether prescribed or prognosti- 
cally computed, that realistically capture inter-annual LAI 
variability. In particular, two characteristics of the LAI 
seasonal cycle affect the extent of LAI control: variability in 
maximum LAI and in the timing of spring green-up and 
autumn senescence. Vegetation types [e.g., DeFries et al., 
1995] that exhibit significant inter-annual variability in both 
green-up time and maximum growing season LAI are espe- 
cially important to represent accurately. Efforts to assimilate 
remotely sensed phenology data [e.g., Stockli et al., 2008a, 
2011] should concentrate on improving phenological esti- 
mates of vegetation with these characteristics. 

[ 49 ] While the interplay between LAI and carbon fluxes is 
particularly important in relation to global carbon emissions, 
the connection between LAI and water fluxes is more rele- 
vant to the simulation of regional climate in GCMs. Pitman 
et al. [2012], in a review of processes needed to improve 
regional climate representation, identified crop phenology 
as one potentially important factor. Although our findings 
show that the controls of LAI phenology are not as important 
for evapotranspiration as they are for transpiration, researchers 
have highlighted the importance of realistic partitioning of 
evapotranspiration into transpiration, soil evaporation, and 
canopy evaporation within GCMs [e.g., Lawrence et al., 
2007]. Importantly, we find that at sites with the greatest 
variability in LAI (S and W sites), phenology has a significant 
control on evapotranspiration timing. Such timing is important 
as it can affect the seasonality of model-generated precipita- 
tion and other atmospheric variables. 

[ 50 ] We find that the correlations between LAI and atmo- 
spheric variables do impact the extremes of carbon and water 
fluxes, especially at two of our sites (NE and FW sites). This 
finding suggests that the representation of LAI correlations 
with atmospheric variables is needed, highlighting the poten- 
tial value of dynamical vegetation modeling. This is particu- 
larly relevant given that climate models are increasingly 
needed for investigations of, for example, temperature and 
precipitation extremes in a changing climate. 

[ 51 ] Our results must be qualified by several limitations of 
our experimental design. The sensitivities shown are for a 
specific land model and would presumably be somewhat 
different for other land models; in this sense, we emphasize 
our first-order findings rather than quantitative details. Also, 
only 17 years of data were available for the statistical analy- 
ses. These limitations present opportunities for additional 
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research with, for example, a recently released 30 year LAI 
data set (http://cliveg.bu.edu/modismisr/index.html) used in 
a coordinated experiment with several of the latest dynami- 
cal models. 
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